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We investigate the time-dependent, coherent, and dissipative dynamics of bound particles in sin- 
gle multilevel quantum dots in the presence of sequential tunnelling transport. We focus on the 
nonequilibrium regime where several channels are available for transport. Through a fully micro- 
scopic and non-Markovian density matrix formalism we investigate transport-induced decoherence 
and relaxation of the system. We validate our methodology by also investigating the Markov limit 
on our model. We confirm that not only does this limit neglect the coherent oscillations between 
system states as expected, but also the rate at which the steady state is reached under this limit 
significantly differs from the non-Markovian results. By a systematic analysis of the decay constants 
and frequencies of coherent oscillations for the off-diagonal elements of the reduced density matrix 
under various realistic tunneling rate anisotropics and energy configurations, we outline a criteria 
for extended decoherence times. 
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I. INTRODUCTION 

It is a fundamental feature of quantum theory that the 
dynamics of an isolated system will follow a unitary evo- 
lution, and thus be fully reversible. In practice, most 
quantum systems are influenced by uncontrolled and in- 
evitable interactions with an often incoherent environ- 
ment, and this influence will typically destroy this de- 
terministic evolution and result in the rapid vanishing of 
any quantum coherence within the open systemi^ 

A large body of theoretical work has been carried out 
throughout the last half century which has pushed for- 
ward the investigation of the dynamics of open quan- 
tum systems (OQS). The seminal work by Nakajima^ 
and Zwanzig"^ on projection operator methods, as well as 
the works by KrauSf4 Lindbladj^ and Gorini, Kossakowki, 
and SudarshanS, on quantum dynamical semigroups, and 
by Davies? on Markovian master equations, have rigor- 
ously established the conditions for positivity of physi- 
cal probabilities in open quantum systems. Since then, 
Markovian approaches have been successfully used^ to 
study steady state phenomena on open systems where 
the past memory of the system is neglected, or on closed 
unitary systems, where the dynamics are deterministic 
and thus reversible. 

The advent of technological applications for quan- 
tum coherence, such as in quantum information pro- 
cessing and cryptography, as well as the use of ultra- 
fast laser pulse excitations, are requiring resolutions of 
the quantum dynamics of a system down to a fem- 
tosecond timescalejS This transient regime, which carries 
the coherence and relaxation dynamics, cannot in gen- 
eral be described by a coarse-grained Markovian limit, 
and a non-Markovian approach is usually necessary. On 
this point there has been significant progress in recent 
years; the application of non-Markovian theories to for- 
mal physical models in the area of low-dimensional dy- 
namical quantum systems has been extraordinarily fruit- 



ful. For example, electron transport by means of full 
counting statisticsriii^ current and shot noise analysis 
by diagrammatic techniques)^ qubit dynamics in the 
presence of 1/f noiseji^ decoherence of qubit systems by 
means of effective Hamiltonians for the reservoir elec- 
tron spin dynamics in quantum dots,!^ and decoherence 
in ballistic nanostructures by projector operator tech- 
niques)^ and transient currents through quantum dots 
by means of the reduced density operator formalism ^iS, 

Formalisms based on non-equilibrium Green's func- 
tions (NEGF) have also been successful in the analysis 
of much phenomena in mesoscopic systems J^^— Notably, 
work in both elastic and inelastic transport in quantum 
dots has been carried out by means of NEGF or deriva- 
tives of this formalismi22i2^ However, NEGF is inherently 
a closed-system formalism where the system has Hamil- 
tonian dynamics,^* and thus does not account for irre- 
versibility arising from interactions with an unseen, un- 
known, or otherwise intractable environment. Promising 
advances have been made to extend NEGF to OQS, for 
example, by treating the environment as a correction to 
the system's self-energy)^ or by calculating two-time cor- 
relation functions,^- effectively separating the time scales 
into transient and steady-state regimes. These and other 
approaches constitute significant progress toward a full 
non-Markovian open system analysis within the NEGF 
formalism. ' 

Even so, when the evolution of the many-body sys- 
tem states themselves are sought, one must calculate the 
Green's functions and then extract the density matrix, 
which in the non-equilibrium case may not be uniquei^ 
An alternative formulation, one we adopt in the present 



t For a rigorous compendium of the theory, developments, and 
general methods to approach OQS, see the excellent presenta- 
tions by Alicki and Lendi^ Gardiner and ZoUer^i^ and Breuer 
and Pettrucionai. 



2 



paper, is to work with the reduced density matrix (RDM) 
directly. This arguably provides a more intuitive descrip- 
tion of the dynamics of non-Markovian OQS, and it read- 
ily yields the actual states of the system and quantum 
coherence between them. 

Although the works cited above have covered a tremen- 
dous amount of ground, to our knowledge there has not 
been a systematic analysis of the effects of diverse en- 
ergy and coupling regimes on the quantum coherence in 
an OQS under the presence of transport-induced relax- 
ation. Thus, the aim of this paper is to investigate the 
effects of varying energy and coupling parameters on the 
non-Markovian dynamics of a single few-electron open 
quantum dot under the RDM formalism. We focus on 
the sequential transport regime where two or more chan- 
nels are available for transport, and solve the RDM for 
the system for the occupation probabilities and quan- 
tum coherence to determine energy and coupling regimes 
leading to larger decoherence times. 

This paper is outlined as follows: In section|lI]we derive 
the evolution equations for the density matrix of the OQS 
in the Born approximation, and a brief description of 
the Markov limit is given for completeness. Section IIIII 
develops the non-Markovian transport theory for a single 
quantum dot nanostructure coupled to semiinfinite leads. 
The tunnelling Hamiltonian for the transport model and 
the assumptions made for the system and reservoirs are 
described in Sec. lIII Al Section lTlI BI presents in detail the 
derivation of the analytical expressions for the transition 
matrix (the memory kernel) containing the dynamics for 
this model. In Sec. IIVI we illustrate the non-Markovian 
approach by considering the case of a multilevel dot in 
a regime where two channels are available for transport. 
The analytical expressions for the elements of the RDM 
are shown in Sec. IIV Al and results are compared with 
the Markov limit under a variation coupling anisotropy to 
different orbitals in the quantum dot. The off-diagonal 
elements representing the coherence in the system are 
presented in Sec. IIVBI as well as a thorough analysis of 
the effects of varying the dominant energy and coupling 
parameters for this model. A summary of the results of 
the variations of the energy and coupling parameters is 
given in Sec. IIV CI Finally, conclusions are presented in 
Sec. El 



II. DYNAMICS OF THE OPEN QUANTUM 
SYSTEM 

We concentrate on the dynamics of an OQS which does 
not in general follow a unitary evolution, and where an 
inherent irreversibility may arise due to coupling with an 
essentially Markovian (memory less) environment. When 
information is continually exchanged between an OQS 
and an environment whose degrees of freedom are ei- 
ther unknown, intractable, or uninteresting, the (possibly 
mixed) states of the system alone are described not by 
state vectors, but rather by the RDM of the system^^. 



The RDM can be obtained by taking a partial trace over 
the degrees of freedom of the environment on the total 
density matrix for the OQS plus environment, such that 

pit) = Trenv[fT(t)] =Y,Wn{t)\M)){M)\ (1) 
n 

is the RDM of the system, (j{t) is the total density matrix, 
and Wn{t) is the probability that the system state \^n{t)) 
is occupied at time t. 

For Eq. ([T]) to represent a physical system, where the 
diagonal elements describe the occupation probabilities 
and the off-diagonal elements describe the coherence be- 
tween the OQS states, we require pii{t) > 0, as well as 
\Pi,j{i)\ < 1 for all i and j.^'^ 

The partial trace over the degrees of freedom of the 
environment has the benefit of accounting for the influ- 
ence of the environment on the OQS, but it also has 
the drawback of only allowing for limited or average en- 
vironment descriptions. Although partial-trace-free ap- 
proaches have been developed^i, to our knowledge they 
have only addressed the steady state regime where mem- 
ory effects have been neglected. 

Considering an environment with a much greater num- 
ber of degrees of freedom relative to those of the system, 
and for weak OQS-environment coupling in a sequential 
tunneling regime (Born approximation), we can assume 
that the system has negligible effect (no back-action) on 
the environment, that the system is not correlated with 
the environment, and that the environment is essentially 
in equilibrium at all times <^ In this case, the total den- 
sity matrix for the OQS plus environment can be written 
as 

a{t)^ p{t)®Q{Q). (2) 

where q{Q) is the environment's equilibrium density ma- 
trix. 

The dynamics of the RDM in the Born approxima- 
tion can be obtained by an iterative integration of the 
Liouville-von Neumann equation to obtain the time de- 
pendence for higher order processesji 

p{t) = - C{t)p{0) + f dt'C{t)C{t')p{t'), 
Jo 

(3) 

where C{t)0 — iti~^[H{t), O] is a Liouville superoperator 
for any operator O and Hamiltonian H{t). 

The matrix elements of the above generalized master 
equation for the RDM in the interaction picture can be 
written in the forrn^ 

Pabit) = V / dt'pUt')Rabcd{t - t')e'^^^^^'' , (4) 

where, Rabcd{t — t') is the memory kernel characteriz- 
ing the system, the environment, and their coupling, and 
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where ')abcd = t^ab-tOcd, with ujab = {Ea-Eb)/h denoting 
energy differences between OQS states. 

The system of coupled integrodifferential convolution 
equations, Eq. is solved by transforming it into an 
algebraic system of equations in Laplace space, 



E 

cd 



yVabcdis) Pab{s) = PcdiO). 



(5) 



where, 



Wabcd{s) = ScaSbd (s + iuJab) " Rabcdis) 



Pab{t) = Pab{t)e 



Rabcd{t) = Rabcd{t)e 



For a given set of basis states, an inverse Laplace trans- 
form of Eqn. ([S]) yields the time dependent solutions in 
the Heisenberg picture. 

Analytic solutions in Laplace-space can be obtained 
for only a few OQS states (~ 2 — 4) , since the computa- 
tional effort rapidly increases with the number of avail- 
able transport channels. For larger numbers of channels, 
we adopt a numerical approach to solve the linear system 
in Eq. ©. 

In the regime of sequential tunnelling transport 
through a quantum dot containing uncorrelated elec- 
trons, the important transport channels are those single- 
particle system states whose energy lies within the 
bias window. These single-particle states define the pos- 
sible dynamical many-body states of the system as those 
involved when the single particle states are empty or oc- 
cupied. Thus, for k channels, a minimum of 2^ many- 
body system states are required (for empty or occupied) . 
Furthermore, (2'^')^ density-matrix elements are required 
to describe the population probabilities as well as the co- 
herence between the states. The temporal evolution of 
the density-matrix elements themselves are governed by 
{2^)^ transition tensor elements. In practice, symmetry 
considerations can reduce this number only by a constant 
amount (for details see Ref. ^ pg. 292). 

For the results in the following sections we perform the 
inverse Laplace transform in the form of a Bromwich in- 
tegral^ over a contour chosen such that all singularities 
are to the left of the contour line. The integration is 
numerically performed'^'^ at each time step using a com- 
bination of a 25-point Clenshaw-Curtis and a 15-point 
Gauss-Kronrod integration method. 



Markov Limit 

We validate the use of a non-Markovian approach by a 
brief comparison with the long-time Markov limit This 
approximation assumes that the environment correlation 
functions vanish at such a fast rate that the reversibility, 
and thus memory of the OQS is essentially destroyed.- 
In such a case, p{t) becomes local in time. 



Pab{t) -^Y.P^d^'^'^ / dt'Rabcd{t-t')e'^' 



c.d 



cdt 



(6) 
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FIG. 1: (color online) (a) Schematic representation of the 
single QD model coupled to source and drain reservoirs. 
Two transport channels are present within the bias win- 
dow, (b) Representation of a orbital asymmetry where edge 
(p-type) orbitals are more strongly coupled to the reservoir 
compared to core (s-type) orbitals. 



For time intervals t ^ t' much greater than the environ- 
ment's correlation time r, the correlation functions for 
environment operators rapidly become uncorrelated and 
decay to zero, {Flit - t')Fk') « {Fl{t - t')){Fk') ~ 0. 
Since the environment is uncorrelated beyond this t—t' 3> 
T time, there are no contributions to the dynamics, and 
the upper integration limit in ^ may be extended to in- 
finity with negligible error in the calculations. The lower 
integration limit of Eq. (|6]) can be taken to negative in- 
finity, and the memory kernel becomes a delta function 
of the channel energies. This leads to a Fermi's Golden 
rule for transitions with strict energy conservation. Thus, 
the coupled set of integro-differential evolution equations 
become a coupled set of first order ordinary differential 
equations. 



Pab{t) = '^Wabcd Pcd{t), 

c,d 



(7) 



where Wabcd is a time-local transition rate. 

The Markov approximation is typically valid only in 
the long-time limit and leads to a zero-average of the off- 
diagonal elements of the OQS, thus the transition rates 
Wabcd can be written as Waabb = Wab, since only transi- 
tions between diagonal elements become relevant. 

In the following sections we verify the validity of the 
Markov limit only once the system has established a 
steady stated, and that both the approach and relaxation 
times to the steady state in this limit are unreliable. 



III. DYNAMICS IN A SINGLE QUANTUM DOT 
A. Transport Model 

We consider a sequential transport model for a single 
quantum dot coupled to semiinfinite leads as depicted in 
Fig. [1^. The total Hamiltonian describing this coupled 
system plus environment is given by. 



H^Hs 



H, 



QD 



Hr 



■Ht, 



(8a) 
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where Hs and Hd are the Hamiltonians describing the 
source and drain leads, respectively. These are taken to 
be non-interacting Fermion systems shifted by the bias 
voltage Vb- 

^siD) = Yj'^'^^i^) ^ \(^^B)dlf^d)(^s{d), (8b) 

s{d) 

with a creation operator for the source (drain), and 
ds(^d) an annihilation operator. 

The Hamiltonian for the quantum dot in Ea. (|8ap is 
given by, 

Hqd = ^(?ICJ,; + eVg)clc^ + Vint, (8c) 

i 

where c| and Ci are system creation and annihilation op- 
erators, the single-particle energies fwji are all shifted by 
the applied gate voltage, Vg, and l^nt is the interaction 
among the confined particles. In the present work we 
set T^nt = 0. Although we expect our qualitative results 
to remain unchanged, the implications of a full coulomb 
interaction in this transport theory may have significant 
consequences and will be the subject of future work. 

The tunneling Hamiltonian describing the coupling be- 
tween the reservoirs and the dot is given by^^i 

'^T^ J2 {n4ck + h.c), (8d) 

k,r—{s,d) 

where T^^''-' is a tunneling coefficient for a particle tun- 
neling from the single-particle state \k) in the dot to the 
source (drain) reservoir. For the present model, we as- 
sume a small level broadening F, and low temperature 
such that F <C ksT <^ AE. In this case the reservoirs 
are in their respective non-degenerate ground states. 



B. Memory Kernel 

In order to obtain the system dynamics for this trans- 
port model, we derive the memory kernel Rabcd{t) ap- 
pearing in Eq. By evaluating the Liouville super- 
operator appearing in Eq. ([3]) in the interaction picture, 
and within the weak interaction limit to second order in 
Eq. ((8d1) . the time evolution of the RDM becomes, 

p{t) ^ - '-TrB[nT(t),p{Q)gm 

-■^^ dt'TrB[nT{t),[nT{t'),p{t')e(m- (9) 

The tunneling Hamiltonian (j8d|) can be rewritten in 
the form, 

Urit) = [Fumlit) + Flmkit)) , (10) 

k 



where Qi{t) = e'-^Qo*Cie-'-^Qo*, and where 

F,{t) = ^ ri;,i,(i) + J2 T^,dLdit) (11) 

s d 

is a generalized reservoir operator that takes into account 
both source and drain reservoir states, where Ls(^j^{t) = 

Noting that the interaction Hamiltonian has zero ex- 
pectation with respect to the equilibrium ensemble of the 
environment, we have that 

TrR[Fk{t)g{0)] = {Fk{t))=0, (12) 
{F,{t)F,,{t'))^{Fl{t)Fl,{t'))^0. 

Thus, the odd moments of Ht as defined in Eq. (fTTj) 
with respect to the density matrix for the equilibrium 
reservoir will all vanish. 

With Eq. (in]), the evolution equation (Eq. ©) be- 
comes 

{{Q.m]it')pit') - Q]it')pit')Q:it)) {FHt)F,{t')} 

+ {p{t')Q^m]{t') - Q^{t)p{t')Q]{t')) {F,{t)F]it')) 

+ {pit')Q,{t')Qlit) - Qlit)pit')Q,it')) {F^imit)) 

+ {Ql{t)Q,{t')p{t')-Q,{t')p{t')Ql{t)) {F,{t)F]{t'))} , 

(13) 

A matrix element of the first term on the right of (fT3)) 
(where there are 8 terms in total) can be rewritten in the 
following way, 

{m'\Q,{t)Ql{t')p{t')\m) = 

= J2 Sn^n'{m'\Q,\k){k\Ql{T)\n){n\p{t')\n')e''''^>, 

(14) 

where = [E,,,, - E,, + -B„' - E„^), r = t - t' , and 

where the RDM elements have been rewritten as, 

n 

Terms 3,5,7 in the right side of (1131) can be written in 
a similar way. 

The elements of the second term in (1131) can be rewrit- 
ten as, 

{m'\Q,{t)p{t')Ql{t')\m) = 
J2 ^rn^"' i^'\Qi\k) {k\Ql ir)\n) {n\p{t')\nY'''^> , 

n,n' ,k 

(16) 
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and similarly for terms 2,6,8 in p3)) . 

The reservoir operator Fi has no explicit time de- 
pendence; we can therefore write, {Fl{t)Fp{t')) = 
(FlFpir)), and (^^(iOi^aW) = (F^Mi^.). 

Finally, for a given quantum dot state, we allow cou- 
pling with all available reservoir states, while treating the 
reservoir energy as a continuum. That is. 



where denotes the reservoir energies, Nr and /Lt^ are 
the 2D density of states and chemical potential, respec- 
tively, for reservoir R, and (/)t and (f>B are the top and 
bottom energies of the band respectively. With the re- 
definitions (Il4l)-([ni), we rewrite Eq. (US]) as. 




(17) 



t'A'n{r){F,F}{T)) + Qi,^Q'l,Jr){F] F,{t)) 



J2 \Qt Q\^,, (r) {Fl {t)F, ) + QiaUir) {F. {r)Fj 



E^™'" \Q'^'kQlrn{r){F,F^{r))+Qi,,Q'^^ir){F]F,{T)) 



(18) 



where, Qj,,fc(r) = {n'\Q,{t ~ t')\k). 

Equation (jl8p is the generalized master equation for 
the system, and describes the time evolution of the RDM. 

By making use of Eq. PT|) and PT|) . we derive the 
reservoir correlation functions appearing in (jl8p to be 

{Fl{t)Fp) = K^pn'^^^.it); {FpFiit)) = Klpn%^^.{t) 

{FlF^{t)) = Kiji%;^.{t)- {F^{t)Fl) = Ki^n%;^.{t), 

(19) 

where we have defined 

^-^{t) ^ ^ — ^ — , ^ ^ {t:tp\ 

representing the dynamics due to the interaction of the 
environment with the quantum dot, where uJxa — {Ex ~ 
Ea)/h is a frequency, {Ta)^ is the tunneling coefficient 
appearing in Eq. (jSd]) . Nr is the density of states for 
the 2DEG environment arising from the continuum of 
reservoir states, and t is time. 

Comparing Eq. ([18]) with Eq. (g]) and using Eq. (fT9|) 
we finally arrive at the microscopically derived expression 
for the memory kernel for the present model 

aJ3,r 

(20) 



where 

Khcd = 'Ac\\c){d\cfi\h) ~&aMcac\\d) 

represents the allowed system transitions for the present 
model, with the indices a, &, c, and d denoting many- 
body states; a, and /3 denote single-particle states, and 
r G {source^ drain}. For the present analysis we assume 
a large band limit where the top {4>t) and bottom ((/)_b) of 
the band are taken to be effectively positive and negative 
infinity respectively. 

The dynamics of the density matrix can be calculated 
for a given set of basis states, where Eq. ([^0]) allows 
for independent tuning and analysis of tunnelling rate 
anisotropics due to asymmetries in the source and drain 
barrier widths, and due also to asymmetries in the cou- 
pling between orbitals with differing angular momentum. 

In the following sections we treat the case where two 
transport channels are available, and analyze the evolu- 
tion of the system states when the energy and interaction 
parameters are varied. 



IV. THE 2-CHANNEL SYSTEM 

We present results for both the diagonal and off- 
diagonal elements of the RDM, under the sequential 
transport model and formalism outlined in the preced- 
ing sections. For definiteness, we consider a quantum dot 
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|0) |1) |2) |3) 

FIG. 2: (color online) Schematic representation of the sys- 
tem under consideration. Four system states are considered 
depending on the relative occupation of the single-particle 
orbitals \a) and |/3), with energies Ea and respectively. 
Each orbital can in principle have different tunnel couplings 
to the leads. 



volves a minimum of 2*^ states. We denote by |0) the 
{N — l)-particle ground state of the system, jl) and |2) 
denote ground-states of the N and {N + l)-particle sys- 
tem respectively, and |3) denotes the first excited state 
of the A^-particle system. The four basis states of the 
system are thus defined as, 

|0) = |(iV-l),.,.), \l) = \iN),.s.), 

|2) = |(7V + 1),.,.), |3) = |(7V)e...). (21) 

A schematic representation of the system is shown in 

Fig.m 



with N confined particles, and with two tunnelling trans- 
port channels available within the bias window. Each 
channel involves a fluctuation in the particle-number be- 
tween N and A''±l, and we choose them to involve either 
the ground or first excited state of the iV-particle system. 
This is an experimentally relevant regime. (See, for ex- 
ample refs. [35ll36| ) 

In general, the availability of k transport channels in- 



A. Population Probabilities 

The occupation probabilities for the system to be in a 
given basis state are obtained from the diagonal elements 
of the RDM. We thus derive the algebraic system Eq. ([5]) 
for the 2-channel system with states defined by Eq. (|21l) . 
The diagonal elements in Laplace space are given by. 



Pnn{s) 



E 



1 r^i r^3 



1 r^ir^j r^k 



Gi r^j r^k /^k r^j 
j^k^k ~ ^j^k^; 



1 - iqc^ 



(~<i (~<3 (~<k 
'-'j'-'k'-'i 



1 r^i r^3 r^kr^l 



(22) 



where S^fe; is the symmetric permutation tensor given and all indices run over the many-body system states 

by, 



^ijkl 



1, ijkl permutation of 0123 
0, otherwise 



defined in Eq. i^. 

(23) In Eq. ([22]) the terms G] and G} are defined as 
I 



G] = D- 



s + iuji^3 - i?i,3, 1,3(5 + S + iuJs^i - i?3,l,3,l(s + i^^3s) 



with, 



G)^D-' 



^l,3,i,i(s + «Wi,3)i?i,i,i,3(s) R3,l,j,j{s + i0J3,l)Ri,i,3,l{s) 



s + ia;i,3 - i?i,3,i,3(s + ioJi,3) s + iwsj - i?3a,3a(s + iw3,i) 



D ^ s ~ Ri^i^i^i{s) 



^l,3,j,i(s + iWi,3)Ei,i4,3(s) R3,l.,i,i{s + i^3,l)Ri,i,3,l{s) 



3 + iuJi,3 - i?l,3,l,3(s + lWi,3) 3 + IW3,1 - i?3,l,3,l(s + iuJ3,l) 
I 



In the Markov limit (Eq. ^) the system of equations domain to obtain, 
for the two-channel case is analytically solved in the time ^ ^ ^ 

P{t)nn = — 



'^ijkl Pijkl 



(25) 
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where, P^jki = WijWjkWki'E.ijki, with E^jki given by 
Eq. ([25)1 . and the transition rate Wij given by, 

^ J ~ OO 

X [S (e^; - Wj-,) 8 (/i'" - e^) + (5 (wy - e^) 8 (e^ - /i'")] 

where r — source, drain, S is the Dirac-delta function, 
and 8 is the Heaviside function. 

As a case study to vahdate the non-Markovian ap- 
proach, we analyze the effects of an orbital anisotropy 
on the dynamics of the system; we define edge and core 
system orbitals by the strength of their coupling to the 
reservoirs owing to the detailed shape of the wave func- 
tion^i In 2-dimensional parabolic confinement, for exam- 
ple, the particle's wave function is composed of associated 
Laguerre polynomials where the core orbitals are weakly 
coupled to the leads owing to a poor (s-type) overlap 
with the reservoir states, whereas edge orbitals (p-type) 
are more strongly coupled to the reservoirs since these 
wavefunctions penetrate more deeply into the confining 
electrostatic barriers (see illustration Fig.lTJj). 

In relation to Fig. [21 we denote orbital |a) an edge 
orbital and |/3) a core orbital. Further denoting the or- 
bital transmission by Ta and Tp , we can define an orbital 
anisotropy parameter as e = 1 — Tp/Ta, with < e < 1. 

The following energy parameters are used in the sub- 
sequent analysis: applied bias Vbias = 6 meV symmetric 
about the Fermi energy epcrmi = 30 meV, with two chan- 
nels within the bias window, Ea ~ eFcrmi + 1 mcV and 
Ep = epcrmi - 1 nieV. 

The non-Markovian system of Eq. (P^ . is evaluated 
numerically by performing an inverse Laplace transform 
as described in Sec. HIl thus obtaining results for the di- 
agonal elements of the RDM in the time domain. 

The time evolution of the diagonal elements of the 
RDM in the large bandwidth limit are shown for increas- 
ing orbital anisotropy in Fig. [3] The set of figures on 
the left (Fig. [3l(a)-(d)) are the results obtained from the 
Markov limit by solving Eq. ([7]) , while the figures on the 
right (Fig. [31(e)- (h)) present the non-Markovian results 
based on Eq. Four features in particular are evident 
in the non-Markovian results which we now discuss. 

First, in the non-Markovian results, we observe high- 
frequency oscillations representing the coherent evolution 
of the occupancy between states with the same particle 
number, |1) and |3) (represented by RDM elements pii(t) 
and P33(t) see Eq. (I^ l. As expected these oscillations 
are not present in the Markov limit. 

Second, for small couplings to the core orbitals {Tp ^ 
Ta) the four probabilities couple into two distinct pairs, 
states |1) and |2) — the N and {N + l)-particle ground 
states — are strongly coupled, as are states |0) and |3) — 
the (A^ — l)-particle ground state and the A''-particle first- 
excited state. We can understand why this occurs with 
recourse to Fig. [51 The transition between states |1) and 
|2) is through an edge orbital (with energy Ea) and this 
coupling is stronger than the transition between |1) and 



(a) PooW 

E=l Pll<' 

P22(t) 

I P33(') 

; 

J- 


(c) PooW 

: ^ = l PiiS 

! P22W 
\ P33W 







(b) PooW 

r E = 0.16 PllW 

1 P22(') 
[ P33(') 


1 (f) PooW 

E = 0.16 PiiW 

1 P22W 

P33W 

1 


(c) PooW 

E=o.o6 Pii<;! — - 

1 P22W 

P33W 


1 (g) PooW 

E=0.06 PnW 

P33W 




1 1 1 1 

(d) PooW 

E = 0.04 PllW 
P22W 

^ P33W 


1 (h) PooW 

E=0.04 PllW 

\ P22W 

P33W 


1 1 _] 1 





01234012345 
time (picoseconds) time (picoseconds) 

FIG. 3: (color online) Markovian and non-Markovian time 
evolution of population probabilities in a quantum dot with 2 
transport channels and four states. (See Fig. ([2]).) The plots 
are for symmetric source and drain tunnel barriers, and vary- 
ing orbital asymmetry. We assume a 6 meV bias symmetric 
about a 30 meV Fermi energy, and two transport channels 
at energies epormi ± 1 meV. Plots (a) through (d) are results 
for the Markov limit (Eq. ^), whereas plots (e) through (h) 
present results for the Non-Markovian theory. 



|0) which involves a core orbital (with energy Ep). Sim- 
ilar arguments apply for the coupling between states |3) 
and |0) (strong coupling) and between |3) and |2) (weak 
coupling). This pairing of probability curves is also ob- 
served in the Markov limit, but again without oscillatory 
behaviour. 

Third, in the steady state (t — >■ 00), all occupation 
probabilities tend to the same value of 1/4. This is seen 
regardless of the tunnelling strengths of core and edge 
states as long as these rates are non-zero. The equal 
probability of 1/4 for each level can be understood as 
being due to the symmetric barriers between the dot and 
the source on the one hand, and between the dot and the 
drain on the other. Since these barriers are symmetric, 
any level will have an equal probability of being either oc- 
cupied or unoccupied. Therefore in the steady state and 
under symmetric barriers the levels will be, on average, 
occupied and empty for the same amount of time. This 
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result shows agreement with the Markov hmit, where as 
expected all dynamical states would share equal proba- 
bility in the steady state under symmetric barriers. 

Fourth, although we observe that all four states are 
equally occupied in the steady state, the time to actu- 
ally reach the steady states does depend on the relative 
couplings to the core and edge orbitals. As expected, we 
observe a disagreement with the Markov limit, where the 
Markovian results show that although the core or edge 
states probabilities couple at a much shorter time, the 
overall steady state is actually reached much later that 
in the non-Markovian results. 

The fact that the time taken to reach the steady state 
increases as the tunnelling to the core state decreases can 
be understood as the effect of decreasing the available 
tunnelling pathways. Fewer pathways available means 
that the system takes longer to reach the steady state. 
This behaviour is exaggerated in the Markov approxima- 
tion. In the limit of zero tunnelling to the core orbital, for 
example, states |3) and |0) will never become occupied, 
and the remaining two levels will each reach an occupa- 
tion probability of 1/2 rather than 1/4, as expected by 
preservation of the trace. In this case the Markov limit 
agrees with the non-Markovian results. 

Although much of the recent work in the area of 
non-equilibrium quantum dot systems has made use of 
Markovian-type approximations,—"— it is fundamental 
to acknowledge the strict long-time limit requirement of 
these types of approximations in general, and thus their 
unreliability on the transient behavior of the system^iS 
especially when considering coherence properties of the 
system at short time scales. The disregard for the mem- 
ory of the system in effect neglects most of the details 
of the coherent transient dynamics, and their use should 
only be invoked for coarse grained phenomena. Nonethe- 
less, as expected, the Markov limit does yield reasonable 
results at sufficiently long times ^i^S when the system has 
reached a steady state, or when considering intermediate 
time scales of averaged behavior, such as the average to- 
tal current through a system. 

Focusing now on the transient results shown in 
Fig. [3lje-h), we have extracted the positions and mag- 
nitudes of the peaks and fit^ these to a decaying expo- 
nential of the form Ae~^*. Results for the decay con- 
stant 7 for each value of £ = 1 — Tcoro/Tcdgc are shown in 
Fig. m^a). We observe a slower decay rate as the orbital 
asymmetry (e) is increased. We interpret this result by 
considering the change in relative occupation probabil- 
ity of the core channels relative to the edge channel as 
tunnelling to and from them is increasingly suppressed. 
In such a case the relative occupation of core states also 
decreases and any contribution from these states to the 
coherent evolution of the system will delay the overall 
decay to steady state, thus increasing the decoherence 
time. 

In a similar manner we have extracted the frequency of 
oscillation as e is increased, where the results are shown 
in Fig. 21 The linear decrease of the frequency with an 
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FIG. 4: (color online) (a) Variation in the rate 7, of the decay 
of the oscillation of pi,i. (b) Variation in the oscillation fre- 
quency of pi.i{t) as the orbital asymmetry is increased (e = 1 
is maximum asymmetry). 



increase in e can be explained by again consider a de- 
crease in the core channel's occupation probability as e 
is increased. Therefore, a coherent superposition between 
states dependent on both core and edge channels will also 
exhibit the relative decrease in probability from the core- 
orbital state. 

In the following section, we look at the coherence 
elements of the RDM under variations of the orbital 
anisotropy as well as of other relevant coupling and en- 
ergetic parameters. 



B. Coherence 

The evolution of the quantum coherence between sys- 
tem states with the same particle number,— the Hilbert 
coherence, is coupled to the evolution of the occupa- 
tion probabilities through the off-diagonal elements of 
the RDM. For the present 2-channel case, these coher- 
ence elements are derived in Laplace space to be. 



/5?3 



J2n=0 Pnn{s)Rl3nn{s + lOJls) 



s + iujrs - i?i3i3(s + iwia) 



(26) 



where the diagonal elements p„„(s) are given by Eq. (|22p . 
and aJi3 is the energy difference between the A^-particle 
states |1) and |3). 
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We numerically evaluate the Bromwich integral of 
Eq. (pS)) to obtain the real and imaginary parts of the 
off diagonal element pi3{t) = {l\p{t)\3) in time-space. 
We are interested in the effects of the barriers and en- 
ergy configuration of the system on the evolution of the 
coherence between dot states. In this regard, the impor- 
tant parameters considered are: 1) the orbital anisotropy, 
2) the source/drain barrier asymmetry, 3) the thickness 
of the barriers (given inversely by the strength of the cou- 
pling (Ta)r between reservoir r and dot orbital a), 4) the 
energy level spacing AE between the transport channels 
present within the bias, 5) the applied bias voltage Vbias, 
and 6) the shift of energy levels within the bias window. 
The variations of these parameters are discussed below. 



1. Orbital Anisotropy 

In Figure EJ we present results for the same energet- 
ics and variation in the orbital anisotropy as shown for 
the diagonal elements (Fig. [3]). The Rabi-type coherent 
oscillations observed in Sec. IIV Al between the diagonal 
elements pii{t) and P33{t) are also clearly seen in the 
evolution of pisit) as expected. 

Investigating the effects of this orbital anisotropy on 
the evolution of the system, we focus on the decay con- 
stant and frequency of oscillation of the off-diagonal ele- 
ments of the RDM as (e = 1 — Tcorc/Tcdgc) is increased. 
In Fig. EKa) the decay of the probability to a steady state 
is fitted^ by an exponential function with rate 7. We ob- 
serve a slower decay rate as the orbital asymmetry factor 



Re(p„(l)) 

; „„ ,-.0.0 Im(P,j(l)) - 


ReCPisM) 

I'l (b) E = 0.6 Im(Pi3(t)) 






Re(Pu(l)) 

Ic) E = 0.9 Im{pi3(t)) . 


Re(p,,(l)) 

(il) 1 = 0.58 ImCPud)) 







0.5 1 1.5 2 2.5 0.5 i 1.5 2 2.5 3 

time (pico.seconds) time (picoseconds) 



FIG. 5: (color online) Time evolution of the real and imag- 
inary parts of the off-diagonal elements of the RDM, pi3(t), 
representing the Hilbert coherence of the system. The plots 
are for symmetric source and drain tunnel barriers, and vary- 
ing orbital asymmetry. We assume a 6 meV bias symmetric 
about the Fermi energy, and two transport channels at ener- 
gies epormiiS meV. The coupling strength for the edge channel 
is set at ftT^jfcdgc) ~ 0-5meV , while the coupling to the core 
channel is varied. 
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FIG. 6: (color online) (a) Variation in the decay constant, 
7, of the oscillation of pi3(t), as the orbital asymmetry is 
increased (e = 1 is maximum asymmetry), (b) Inset, variation 
in the oscillating frequency of pi3(t). 



(e) is increased. As in Sec. IIV Al we interpret this result 
by considering the change in relative occupation prob- 
ability of states dependent on the core channel relative 
to states dependent on the edge channel, as tunnelling 
to and from the core channel is increasingly suppressed. 
In such a case the relative occupation of core states also 
decreases and any contribution from these states to the 
coherent evolution of the system will delay the overall 
decay to steady state, thus increasing the decoherence 
time. 

The effect of a variation of the orbital anisotropy on 
the frequency of oscillation is shown in Fig. ^h). We 
again observe a linear decrease of the frequency with an 
increase in e. By the same argument as discussed in 
Sec. IIV Al we consider a decrease in the occupation prob- 
ability of the core channel as e is increased. A coherent 
superposition between states dependent on both the core 
and edge channels will also exhibit the relative decrease 
in probability from the core-dependent state. Thus as 
expected, the coherent evolution between states with the 
same particle number is described by the off-diagonal el- 
ement pi^it). 

2. Source/Dram Barrier Asymmetry 

Keeping now the orbital anisotropy fixed and varying 
the asymmetry between the source and drain barriers, we 
observe that as the thickness of one barrier relative to the 
other is arbitrarily changed, there is a always an increase 
in the decay rate 7 (see Fig. [7^). This can be under- 
stood by considering the contribution to the decoherence 
rate due to various asymmetries: In the case where the 
drain barrier is thicker than the source barrier, tunnelling 
through the drain is suppressed relative to that through 
the source barrier. Therefore decoherence due to rapid 
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FIG. 7: (color online) (a) Variation in the decay constant, 
7, of the oscillation of the coherence element pi3(t), as the 
asymmetry between the source and drain barriers is increased 
(Tdiain/Tsource 7^ !)• (b) Inset, Variation in the oscillating 
frequency of Pl3(^). 
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Coupling Strength (meV) 

FIG. 8: (color online) Variation of the decay 7 (a), and 
frequency of oscillation cj (b), as the overall coupling be- 
tween orbitals and reservoirs is increased. The overall cou- 
pling is changed by varying the strength of the tunnelling 
rate Tcdgo and by keeping a fixed orbital anisotropy of e = 

ore /Tcdgc — 0.3 



tunnelling out of the dot is also suppressed. In the op- 
posite case where the source barrier is thicker, tunnelling 
into the dot is suppressed relative to tunnelling out and 
thus, decoherence due to fast injection into the dot is 
also suppressed. Clearly, the decoherence will be great- 
est in the completely symmetric case since the contribu- 
tions from both electron injection and extraction will be 
greatest. 

Similarly, we find that the frequency of oscillation 
of the coherence term {piait)) is also maximal in the 
completely symmetric case, and decreases otherwise (see 
Fig. [71d). This is expected since we can consider the fre- 
quency of oscillation to be proportional to the energy 
splitting due to the coupling with the reservoirs. Note 
that the shape of the curve in both the decay constant 
and frequency plots is quadratic with the relative thick- 
ness of the barrier. This is a consequence of the fact that, 
to leading order, dissipative effects are quadratic in the 
system-bath couplings. In the limit of infinitely thick 
barriers, there is no oscillation since tunnelling though 
the barrier has been completely suppressed. 

3. Coupling Strength 

We now investigate the dependence of the coherence 
on the coupling strength between the reservoirs and the 
available orbitals in the dot while keeping the asymme- 
try between the barriers fixed. For simplicity, we also 
keep the ratio of the core-orbital tunneling rate to the 
edge orbital tunneling rate fixed such that the coupling 
strengths are varied at the same rate. 

Figure [S] presents the change in the decay constant 
and the frequency of oscillation of the coherence element 
Pi3{t) as the coupling strength is varied. We observe 



that both the decay constant and frequency closely fol- 
low a quadratic fit on the coupling. As mentioned in 
the discussion of the results of barrier asymmetry, this is 
in agreement with the coupling effects being quadratic to 
leading order. We expect this to be a reasonable result in 
cases where sequential transport is the dominant process. 



4-. Energy Level Spacing 

Keeping now the barriers fixed, as well as the applied 
bias voltage, we look at a variation in the energy level 
spacing AE = Ea — Ep between the A^-particle states |1) 
and |3). (All previous results considered AE = 2meV.) 
Results for the decay constant and frequency are shown 
in figure [IJa) and (b), respectively. The level spacing 
is varied symmetrically about the Fermi energy. In this 
case we observe a decrease in the rate of decay as the 
level spacing is increased (Fig. This result can be 

understood in the following way: consider an electron in 
the source reservoir which tunnels into one of the avail- 
able levels in the dot, and then tunnels out to the drain 
reservoir. Since the dot states are not eigenstates of the 
full Hamiltonian, the energy of the electron within the 
dot is not well defined and can have a value within the 
energy level spacing, AE. As the electron tunnels from 
the source reservoir into the dot it can have at most an 
energy Smax = eFormi + l/2eybias- For the electron to 
then tunnel into a drain state, it must have an energy 
of at least i?min — epormi — l/2eT43ias + AE. There- 
fore, the available energy states of electrons that can 
tunnel into the dot from the source reservoir is given by 
Emax — E^in = eVhiag — AE. Thus, the greater AE is, the 
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FIG. 9: (color online) (a) Variation in the exponential decay 
constant 7, as the energy level spacing between the |a) and 
single particle levels is changed. Error bars indicate how well 
the decay is approximated as being exponential, (b) Variation 
in the frequency of coherent oscillation oj between many-body 
A''-particle states |1) and |3). 



fewer available energy states for transport into the dot, 
and the longer lived the dot states are. At approximately 
1.5 meV we observe a sudden drop in the decay constant 
leveling again after an increase in AE of 0.5 nieV. We 
investigate this feature by looking at its dependence on 
other parameters, and find that as the coupling strength 
between the reservoir and the dot is decreased, the fea- 
ture decreases (see Fig. [9}d). We also observe that the 
slope of the curve decreases after the drop, indicating a 
slower decrease of the decay constant. 

The dependence of frequency of oscillation of the co- 
herence on the energy level spacing is shown in Fig. |9j3. 
The relevant energy differences giving rise to oscillations 
are the level spacing AE on the one hand, and the level 
energy and the edge of the bias window on the other. 
For small AE, it is the latter which dominates and so uj 
shows weak dependence on AE as shown in Fig. [Hb. For 
larger AE, it is AE itself which dominates and so we see 
UJ growing linearly with AE. 



5. Applied Bias Voltage 

A variation of the applied bias voltage has also been 
investigated, with the results presented in Fig. [101 Here 
Epermi = SOmcV, E^ = 31meV, E/s = 29meV. For a 
symmetric variation of the bias voltage about the Fermi 
energy, we observe an overall increase of the decay con- 
stant as the bias voltage is increased (Fig. [TUk). We 
explain this overall behaviour by considering that as the 
bias is increased, more energy states are available in the 
reservoirs for electrons to tunnel into or out of, elec- 
trons in the dot then have a larger number of tunneling 
pathways, thus increasing relaxations, and decoherence. 




5 6 7 

FIG. 10: (color online) Variation of the decay constant 7 (a,c) 
and of the frequency of oscillation uj (b,d) as the applied bias 
voltage is increased symmetrically about the transport chan- 
nels. In (a) and (b) the transport levels are kept symmetric 
about the center of the bias window, whereas in (c) and (d) 
the transport levels are shifted off-center within the bias win- 
dow by 0.5meV. In all cases, AE — 2meV . 



This overall trend is also seen in the frequency of os- 
cillation of the coherence (Fig. [TOb) where, as the bias 
is increased, the frequency generally decreases. As dis- 
cussed in Sec. IIV B 41 the frequency of oscillation follows 
the magnitude of the level spacing relative to the energy 
difference between the levels and the edge of the bias. 
Therefore, as the bias increases relative to the level spac- 
ing, the frequency of oscillation is expected to decrease. 

Within this overall behaviour, we observe further 
structure in Fig. [TU] in the form of pronounced plateaus 
occurring at specific intervals. This structure can be ex- 
plained by noting that there are low frequency oscillatory 
envelopes (not shown) on the evolution of the RDM el- 
ements given by the interplay of all relevant energies in 
the system. The intervals from one plateau to the next 
are representative of the frequency of these slow envelope 
oscillations. The effect is also seen in the behaviour of 
the high frequency oscillations, where as the bias in in- 
creased, these form peaks with period given by the low 
frequency envelope (Fig. ITUb ) . 

We also find the number and position of these plateaus 
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FIG. 11: (color online) (a) Variation of the decay constant as 
the transport energy levels are shifted within the bias window, 
(b) Variation of the frequency of oscillation as the levels are 
shifted. Here AE = 2meV. 



of these peaks about the Fermi energy, where the overall 
height decreases as the levels increase past the center of 
the window. As discussed in Sec. IIV B 5[ the plateaus 
appearing in the decay coefficient, and the peaks appear- 
ing in the frequency of oscillation are due to envelope 
oscillations of low frequency in the evolution of the RDM 
elements. 

From these results, we can assert that the decoherence 
time is longest when the levels are either centered in the 
bias window, or with at least one level in resonance with 
the chemical potential of one of the reservoirs. 

Similar to the results of Sec. IIV B 5[ we cannot defini- 
tively ascribe a microscopic origin to these plateaus and 
peaks. Analysis of this structure will be the focus of fu- 
ture work. 



C. Summary 



and peaks depends on the position of the levels within 
the bias window. Figures [TUb and fTUH show an increase 
in the number of plateaus as the levels inside the dot 
are shifted off-center by O.lmeV within the bias window 
{Ea = 31.5meV,Ef; = 29.5meV). Therefore the low fre- 
quency oscillations are strongly dependent on the relative 
energy difference between the level spacing and the spac- 
ing between the levels and the edge of the bias window. 

Further analysis of these features will be the basis of 
future work by the authors. 



6. Energy Level Shift Within Bias Window 

As mentioned in the discussion above, the position of 
the transport channels relative to the center of the bias 
window has consequences on the overall coherence be- 
tween system states. As before, we keep the bias fixed 
at 6 meV symmetric about a Fermi energy of 30 meV, 
giving a chemical potential for the source reservoir of 
i'Fermi + ^^ijias = 33 uieV . Similarly the chemical po- 
tential for the drain reservoir is then ii^Fcrmi — jVbias = 
27 meV. By varying both levels simultaneously and keep- 
ing the level spacing at AE = 2 meV, we can then move 
the edge level {Ea) from 29 meV to 33 meV to sweep 
the levels over the entire window. Figure [TT] presents the 
result of this shift of the position of Ea- In a similar 
fashion to the variation of the applied bias, we observe 
structure in both the decay constant and frequency of 
oscillation of the coherence of the system. In Fig. [TTk . 
the decay constant presents two plateaus nearly symmet- 
rical about the Fermi energy with size in the order of 
{0.5) AE. For the frequency of oscillation (Fig. \TTh) we 
see that it presents peaks at about the same energies as 
the edges of the plateaus (similar to the results of varia- 
tion of Vbias in Fig. ITU)) . We also observe an asymmetry 



We now summarize the results of the variation of en- 
ergy and coupling parameters on the decoherence and 
frequency of oscillations of the system. For the sequen- 
tial transport model presented, the following conditions 
may independently increase the decoherence time in the 
system: An increase in the difference of tunnelling rates 
between available channels in the transport window, an 
increase in the barrier asymmetry between source and 
drain barriers, a decrease of the overall coupling strength 
between reservoirs and dot, an increase in the energy level 
spacing between available transport channels, a resonant 
case between available channels and reservoir Fermi en- 
ergies or a relatively small bias voltage, a symmetrical 
positioning of transport levels within the bias window, 
or resonance between a transport level and a reservoir 
Fermi energy. 

Similarly, the frequency of oscillation can be decreased 
by: An increase in the orbital anisotropy between avail- 
able channels in the transport window, an increase in 
the barrier asymmetry between source and drain barri- 
ers, a decrease of the overall coupling strength between 
reservoirs and dot, a decrease in the energy level spacing 
between available transport channels, a resonant case be- 
tween available channels and reservoir Fermi energies, a 
relatively high bias voltage, a symmetrical positioning of 
transport levels within the bias window. There is also the 
case that the frequency decreases in certain intervals as 
the bias voltage is increased or as positions of the energy 
levels are swept over the bias window. Table U presents 
a tabular summary of the above results. 

We thus find that an optimal lifetime for the coher- 
ence between system states may be obtained by a suit- 
able combination of the above energetic and interaction 
regimes. At the very least, by using the above crite- 
ria, it may be possible to confirm the barrier and energy 
parameters in experiment by studying decay rates and 
oscillatory behaviour. 
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Parameter Decoherence-time Frequency 



Level Spacing Increases 

Bias voltage Decreases 

Off-center shift Decreases 

Coupling strength Decreases 

Barrier asymmetry Increases 

Orbital anisotropy Increases 



Increases 
Decreases 
Decreases 
Increases 
Decreases 
Decreases 



TABLE 1: Effects of increasing magnitude of parameters on 
coherence and frequency of oscillation. 



V. CONCLUDING REMARKS 

We have developed a non-Markovian formalism for the 
transient evolution of the reduced density matrix of a 
quantum dot weakly coupled to source and drain reser- 
voirs. Within a tunneling Hamiltonian approach, we have 
analytically derived an explicit memory kernel describing 
the sequential transport dynamics of the system for an 
arbitrary number of transport channels. The results of 
the analysis where compared with the Markov limit for 
a 2-channel (4-state) system, where we verified marked 
differences in the transient dynamics. Apart from the ex- 



pected absence of coherent oscillatory behaviour in the 
Markovian resiilts, it was also found that the rate at 
which the system approaches a steady state differs con- 
siderably between both theories in the regime of highly 
anisotropic tunneling into distinct system orbitals. It was 
found that the Markov approximation shows significantly 
longer time to reach a steady state when the tunnelling 
anisotropy is high, thus confirming its applicability only 
in the long-time limit. Through a comprehensive and sys- 
tematic analysis, the decoherence time and frequency of 
oscillations observed in the non-Markovian results where 
found to depend on both the energy parameters of the 
system as well as the distinct coupling parameters to the 
reservoirs. Finally, distinct regimes have been outlined 
where both of these coherence characteristics could be 
increased. 
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